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ABSTRACT 


THERMAL ANALYSIS AND DESIGN OF A COOLING 
SYSTEM FOR THE NASA AMES MACH 14 WIND TUNNEL NOZZLE 


This report provides the analysis and design of a Mach 14 
converging diverging nozzle wall liner. The analysis 
indicates that: no fin on the coolant side of the nozzle 

wall is optimum, the thermal stresses are dominant, and the 
critical area is very near the throat. The molybdenum alloy 
TZM, with a wall thickness of 2.0 mm in the throat area, 
appears to be the only material capable of meeting design 
requirements. Additionally, cooling water at 2000 psia with 
a flow velocity of 25 m/s in the coolant passages is 
required. 
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2 

A heat transfer area, m 

skin friction coefficient, unitless 
Cp specific heat, kJ/kg*K 

hydraulic diameter, m 
E modulus of elasticity, psi 
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h heat transfer coefficient, W/m *K 

k thermal conductivity of liner material, W/m K 

L length, m 

Ma Mach number, unitless 

Nu Nusselt number, unitless 

p pressure, psi 

P pressure, psi 

Pr Prandtl number, unitless 

q transferred heat energy, W 

r radius, m 

r recovery factor, unitless 
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Re Reynolds number, unitless 

Sy yield strength, psi 

St Stanton number, unitless 

T temperature, °F, °C, Kelvin 

U velocity, m/s 

V velocity, m/s 
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NOMENCLATURE (continued) 

coolant passage width, m 

distance along nozzle wall, m 

distance along nozzle axis of symmetry, m 


SUBSCRIPTS 

stagnation conditions 

air side 

bulk conditions 

coolant side 

inside surface 

outside surface 

straight 

radial component 

recovery 

wall 

axial component 
static conditions 
throat conditions 
circumferential component 
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NOMENCLATURE (continued) 


GREEK 

coefficient of thermal expansion, m/m K 

ratio of specific heats, unitless 

velocity thickness, m 

displacement thickness, m 

potential boundary layer thickness, m 

real boundary layer thickness, m 

temperature thickness, m 

temperature difference, Kelvin 

momentum thickness, m 

Poisson ratio, unitless 

dynamic viscosity, kg/m s 

kinematic viscosity, m /s 

density, kg/m^ 

stress, psi 

energy thickness, m 


CHAPTER 1 


Introduction 

In order to test large aerospace models under hyper- 
sonic conditions the 3.5 foot blowdown type wind tunnel 
shown in Figure 1 was designed in 1957 by NASA Ames engi- 
neers. The supply air was compressed and then heated in a 
gas-fired pebble-bed heater to the necessary stagnation 
conditions. Four interchangeable nozzles were to permit 
airspeeds of Mach 5, 7, 10 and 14. However, unexpected 
problems prevented the Mach 14 nozzle from being used. 

The original Mach 14 nozzle shown in Figure 2 is taken 
from a NASA drawing. The converging section of the nozzle 
employs a water cooling system, while the diverging section 
uses film cooling via injection at the throat. It is 
desired to eliminate the film cooling and switch to a water 
cooling system for the diverging section as well. The 
nozzle is to be operated at an inlet temperature and 
pressure of 3460 °R and 2000 psia, respectively. The 
objective of this research is to design a cooling system 
which will allow the nozzle to operate within its material 
limitations. 


1 





'/■2A/2-OG 7SO- 


OmOtNAL PAGE IS 
or POOR QUALITY 



Figure 2. Cross section of original Mach 14 nozzle 




CHAPTER 2 


Design 

The design of the cooling system for this nozzle began 
with the following requirements: (1) The nozzle is fixed 
and is given in Figure 3. (2) The nozzle profile is sup- 

plied with dry air at 2000 psia and 3460 °R. The next 
requirement was the obvious necessity of having the nozzle 
operate within its material limitations. This encompasses 
two considerations: (1) The material must not exceed its 

oxidation temperature, and (2) the total maximum stress must 
not exceed the yield strength of the material at its 
elevated operating temperature. 

Engineering design is a combination of synthesis and 
analysis. An initial design is analyzed and the results of 
this analysis lead to an improved design, which is again 
analyzed and improved. This process continues until an 
optimum design is attained. 

Past experience showed that the nozzle wall must be 
cooled. The first design thus began with a wall finned on 
the backside as shown in Figure 4. The objective in this 
initial approach was to optimize the fin geometry (dimen- 
sions d^ through d 4 ) . A numerical simulation using a finite 
difference model of the finned wall with appropriate 
expressions for heat transfer coefficients on the air and 
coolant sides, and for a variety of wall materials, revealed 
that no fin was optimal! This conclusion was later 
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In drawing A12286-030 
Dec. 17. 1964 



Figure 3. Mach 14 wind tunnel nozzle geometry 
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confirmed with a mathematical approach involving a critical 
Biot number which is described in detail in Appendix A . 



Figure 4 . initial design of nozzle wall with fins. 

The analysis then focused on a plane wall as shown in 
Figure 5. Preliminary analysis indicated that the most 
critical need for cooling is in a region 5 inches upstream 
to 5 inches downstream of the throat. Thus the design shown 
in Figure 6 was finally arrived at. 



coolant passage 

£ 


nozzle 


wall 


air side 


Figure 5. The plane nozzle wall 
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Note that the liner contains a flanged end downstream 
and a sliding upstream end with a seal. This design allows 
for thermal expansion. The reason for having the downstream 
end fixed and the upstream end free is due to flow con- 
siderations. The free end will present a discontinuity in 
the wall profile. This discontinuity is less detrimental to 
the flow on the converging wall than on the diverging wall. 
The disadvantage associated with this arrangement is that 
the seal of the floating end sees air pressure at about 2000 
psia and nearly 3460 °R, whereas the downstream air pressure 
is only about 4 psia with a temperature well below 500 °R. 
Thus, a special dynamic seal is required in the high 
temperature region. 

Note the taper in the wall thickness from the throat 
section to either end of the nozzle liner. To facilitate 
heat transfer, past experience has indicated that the liner 
wall must be quite thin in the throat region. However, a 
thicker wall is more desirable at the downstream end of the 
nozzle for transition to the flange. Similarly, a thicker 
wall is more desirable at the upstream end where the dynamic 
seal is located. 

Coolant flow is at right angles to the air flow across 
the backside of the nozzle liner. Channel definition is 
accomplished with the bosses or channel dividers shown in 
Figure 6. This design provides predictable heat transfer 
results over the entire cooled portion of the nozzle liner. 




Figure 6 . The proposed Mach 14 nozzle and coolant 
system design. 


coolant exit 
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The coolant passages are well defined, curved rectangular 
ducts 0.4 cm high by approximately 8.0 cm in width. 

In order to maintain a nearly constant coolant velocity 
throughout the cooled portion of the nozzle a series ducting 
system is used as shown in Figure 6. A pressure differen- 
tial between the various coolant passages causes bypass 
coolant flow around the channel dividers thus cooling the 
nozzle wall in this area. High pressure, high flow rate 
coolant is introduced to the nozzle at the downstream sec- 
tion. Any rise in coolant bulk temperature will be minimal 
since the cooling requirement is low in this section. Next 
the coolant is delivered to the throat section where the 
majority of heat transfer takes place. The coolant then 
passes through the upstream flow passage and exits the 
nozzle section. 

With the basic design of the nozzle liner established, 
it is now necessary to optimize the wall thickness profile, 
and to select an appropriate material. Additionally, the 
coolant must be optimized for flow rate, bulk temperature 
and pressure. In the following chapter we develop the 
analyses necessary to choose the best design parameters. 


CHAPTER 3 


Analysis 

The analysis was performed in the following sequence: 

1. Flow analysis. This gives the temperature, pres- 
sure, and Mach number throughout the nozzle. 

2. Thermal analysis. In this analysis the nozzle wall 
temperature distribution is obtained as a function of the 
various design parameters. 

3. Stress analysis. With the known temperature dis- 
tribution thermal stresses are obtained and combined with 
pressure stresses. 

4. Material analysis. This analysis shows the rel- 
ative importance of the physical properties of the candidate 
materials. 

Flow Analysis 

Given the nozzle inlet temperature and pressure the 
local free stream temperature, pressure, and Mach number are 
obtained using compressible flow theory. Given the geometry 
of the nozzle we can solve for the Mach number implicitly 
using the Mach number relation for air ( Y = 1.4) 

A _ 1 (1 +0.2 Ma 2 ) 3 

A* “ Ma 1.728 
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Once the Mach number is known we solve for free stream 
temperature and pressure by 


0 


CO 


1 + 0.2 Ma" 


t 


= P, 


(1 


+ 0.2 Ma 2 )" 3 * 5 


Several values 
are shown in Figure 


of free stream temperature and pressure 

7. 


P 

T 



p * 1056 psla 
T - 2875 °R 


Figure 7. Compressible flow analysis. 


Thermal Analysis 


A cross section of the nozzle wall is given in 


Figure 8. 
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With reference to Figure 8 we can write 

% * % m <3 C • (1 > 

The assumptions are: 

1. steady state conditions. 

2. negligible conduction in the wall parallel to 
its surface. 

3. negligible radiation from the air. 

Using Newton's law of cooling for the first and last 
terms in Equation 1 and Fourier's law of conduction for the 
center term we have 


q 


a 


A h (T 
a a v r 
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w,a 
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( 2 ) 
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A_ h - 
c c 


(T 


w,c 
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(3) 


% * - A r k H 


2^ L k (T Wja - T w-C ) 
ln < r o /r i > 


(4) 


The heat transfer on the air side comes from Bartz [1], 
The complete development is given in Appendix B. Here we 
may write 


h 


a 


< St C P U * >air 


(5) 
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The heat transfer on the coolant side is from Rohsenow 
and Hartnett 12] . The complete development of this analysis 
appears in Appendix C. The result is 


h c = 0.0216 (k/D h ) Re 8 * 8 Pr 1/3 (yu/y^) 0 * 14 . (6) 


We can now obtain expressions for the wall temperature 


on the air and coolant sides (T , T ) by manipulating 

W / cl W f C ** 

Equations 1 through 6. The result is 

k T w,c 


w, a 


T r r i h a + ln(r Q /r i ) 


In (r Q /r i ) + r i h a 


(7) 


and 


w,c 


R T w, a 

in (r Q /r i> T h c r o T b 


r o h c + In (r /~r • ) 


o 1 


( 8 ) 


We now find that 


T = T (h , T ) 
w,a w f a ' a ' w,c 


(7) 


h = h (T ) 
a a w,a 


(5) 


T w,c T w,c ^ T w,a ' ^c^ 


( 8 ) 


h c - h c (T w,c» 


( 6 ) 
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Because of the implicit nature of Equations 7, 5, 8 
and 6 the solution is obtained by the numerical method of 
successive substitution. Details of the computer code which 
performs these operations are given in Appendix D. 

Equations 7 and 8 give the temperature of the nozzle 
wall. For design considerations the maximum wall tempera- 
ture T m _„ must not exceed the corrosion temperature of the 

W f Iucl X 

material. Therefore we define a Margin of Safety for Tem- 
perature 

Margin of Safety _ _ T 

for Temperature w,max corrosion . (9) 


Stress Anali 


The thermal gradient and differential fluid pressure 
across the nozzle liner set up a state of stress. The prin- 
ciple of superposition allows us to calculate these stresses 
independently and then simply add. them together. Thermal 
stress was calculated using the long hollow-cylinder equa- 
tions assuming a logarithmic temperature distribution 
according to Timoshenko [3]. Justification for a loga- 
rithmic temperature distribution as a function of radial 
distance is given in Appendix E. Stresses will be greatest 
at either surface. At the inside radius the circumferential 
and axial components of thermal stress are given by 



( 10 ) 
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and the radial component of thermal stress is given by 



( 11 ) 


At the outside radius circumferential and axial components 
of thermal stress are given by 




o(E 4T 


1 _ 2 r i ln < r n/ r i> ' 


f O - r i 


2 (1 - "X ) ln(r^7r~r 


( 12 ) 


and the radial component of thermal stress is given by 



(13) 


The material properties appearing in Equations 10 and 
12 are dependent on temperature. The pressure stresses are 
calculated using Timoshenko's thick-walled pressure vessel 
theory. It should be noted that thermal stresses are up to 
18 times greater than pressure stresses in the critical 
areas near the throat (see Appendix F) . For this reason the 
axial components of pressure stress, which are a complicated 
function of the nozzle geometry, are omitted. At the inside 
radius the circumferential component of pressure stress is 
given by 





2 P 


c 




t 


( 14 ) 
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where is the local air side pressure and P^, is the 

coolant pressure. The radial component of pressure stress 
is given by 


<r 


r , 1 


= -P 


QO 


• (15) 


At the outside radius the circumferential component of pres- 
sure stress is given by 
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2 

o 


+ 



r 


2 

i 


(16) 


and the radial component of pressure stress is given by 


^ z ,o ~ P c 

The theory of failure used 
and Mitchell [4]. It is called 
of failure. The triaxial stress 


(17) 

here is described by Shigley 
the distortion energy theory 
state is calculated as 


<f y ={l/2 [«T 9 -<f r ) 2 +(^ -Cf z ) 2 +(<f' z + (f r ) 2 ] ] 1/2 (18) 


For design considerations the maximum triaxial state of 
stress must not exceed the yield strength of the material at 
the operating temperature. Therefore we define a Factor of 
Safety for Stress 


Factor of Safety _ 
for Stress ~ <r, 


y ,max 


( 19 ) 
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Before the blowdown process begins, the air side of the 
nozzle wall experiences a hard vacuum, while the coolant 
side is subjected to the full coolant pressure of 2000 psia. 
Failures that are possible during this pre-blowdown condi- 
tion are discussed in Appendix G. 

Material Analysis 

The obvious criteria for material selection are Factor 
of Safety for Stress and Margin of Safety for Temperature. 
There are, however, two additional considerations. The 
first is wall thickness. Preliminary analyses have shown 
that a thin wall is necessary due to thermal stress. But 
from a practical point of view thin walls are difficult to 
fabricate, require delicate handling, and provide little 
allowance for loss of material due to corrosion and erosive 
"sandblasting" that occurs when the wind tunnel is operated. 

The second additional consideration is that the wall 
temperature on the coolant side must not exceed the satura- 
tion temperature of the coolant. Tong [5] has described two 
phase boiling heat transfer and demonstrates that local nuc- 
leate boiling at the heat transfer surface (the coolant side 
wall in this case) actually increases heat transfer up to 
the "boiling crisis" temperature. Boiling crisis occurs 
when local boiling has advanced to the point that a bubble 
layer effectively insulates the heat transfer surface from 
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the coolant causing a radical increase in the surface tem- 
perature. This kind of temperature rise is undesirable and 
may lead to material failure. Using expressions from Tong, 
and with water as the coolant, we find that boiling crisis 
occurs approximately 35 °C above the saturation temperature 
for a large range of coolant pressures. Due to possible 
inaccuracies in the overall computer model temperature cal- 
culations for the current application, boiling is prohibited 
for this design. The upper limit for T w c is therefore the 
critical point for water where the temperature and pressure 
are 374 °C (705 °F) and 22.1 MPA (3206 psia) respectively. 
Saturation temperatures corresponding to lower pressures are 
given in Table 1. 


Table 1 . 

Saturation Temperatures and Pressures for Water 


Pressure 
MPa psig 

Temperature 
°C 0 F 

6.9 

1000 

285 

545 

10.3 

1500 

313 

596 

13.8 

2000 

335 

636 

17.2 

2500 

353 

668 

20.7 

3000 

368 

695 
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Summarizing the criteria for material selection we have: 

1. Factor of Safety for Stress. 

2. Margin of Safety for Temperature. 

3. Wall thickness. 

4. Coolant side wall temperature. 

Preliminary analysis indicates that thermal stress 

dominates over pressure stress and that the thermal stress 
on the air side is slightly greater than on the coolant 
side. The thermal stress on the air side of the nozzle 
liner was given in Equation 10 and is repeated here. 


6 = 


a<E AT 


In U 0 A j) 


( 20 ) 


2 (1 - ) In U 0 /r i ) 


The at is obtained from the heat transfer through the wall. 


2 ITk L (AT) 

ln(r 0 / r i> 


( 21 ) 


Solving for aT 


AT 


q 

2 it k L 


ln < r o /r i> 


Substituting Equation 21 into 20 gives 


<r 


\ 


oc E 

r 

q 

2 r o ln(r o /r i ) 

k 

4 Tr L (1 - Tv ) 

i. 

1 ' r 2 - r 2 

O 1 


( 22 ) 


For a nozzle liner of fixed wall thickness the term in 
the braces of Equation 22 is relatively constant. Although 
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q increases with an increase in k it is only weakly depen- 
dent. We therefore write 

oc E p 

= —IT" 1 



where is a constant. Recalling the Factor of Safety for 
Stress 


F.S . 


S_ 

<r 


So 


F.S . 


Or 


S 


C< E p 

“IT" 


F.S. OC 


S k 
c* E 


(23) 


We conclude from Equation 23 that in order to satisfy 
the Factor of Safety for Stress the material selected must 
have a relatively high yield strength and high thermal con- 
ductivity, and a relatively low coefficient of thermal 
expansion and low modulus of elasticity. Candidate 
materials that possess the above properties at elevated 
temperatures are given in Table 2. 










CHAPTER 4 


Results and Discussion 


The optimum design 

Using the mathematical model developed in the analysis 
chapter we are now ready to select from the following design 
parameters: the material, wall thickness, coolant flow rate, 
coolant pressure and bulk temperature. This will produce a 
final design which gives the optimum combination of Factor 
of Safety for Stress and Margin of Safety for Temperature. 
Additional constraints are a practical wall thickness and 
the wall temperature on the coolant side lower than the 
coolant boiling temperature. Initial examination of the 
design parameters led to the range and combination of para- 
meters given in Table 3. 

The results of the analysis for the combinations of 
design parameters given in Table 3 appear in Figures 9 
through 13. The selected design calls for TZM with minimum 
wall thickness of 2 millimeters, and coolant flow rate of 25 
m/s at 300 K and 2000 psia. The three graphs in Figure 9 
plot various quantities against the nozzle axial location 
for the design case. Figure 9a gives the results of the 
flow analysis; Figure 9b gives the results of the thermal 
analysis; and Figure 9c gives the results of the stress 
analysis . 
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Table 3 


Selection of Design 


Parameters . 



* 


TZM is molybdenum alloyed with titanium and zirconium 
GLIDCOP is copper alloyed with an aluminum oxide. 
AMZIRC is copper alloyed with zirconium. 

WC-103 is a columbium base alloy. 
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The axial locations span from 7 cm upstream of the 
throat to 2 cm downstream of the throat. The reason for 
terminating the information downstream of the 2 cm location 
is two fold: first, the temperatures and pressures are 

rapidly decreasing at this point and beyond since we are 
well past the critical location just upstream of the throat; 
second, the computations become time consuming due the num- 
erical solution of h as described in Appendix C. 

a 

In Figure 9a the pressure, temperature and Mach number 
of the air throughout the nozzle are given. Figure 9b gives 
the heat transfer coefficient, heat flux and wall tempera- 
tures. Notice that all of these parameters peak at about 
0.5 cm upstream of the throat. In Figure 9c the Factor of 
Safety for Stress and Margin of Safety for Temperature reach 
their relative minimums at this 0.5 cm upstream location. 

The effect of varying the nozzle wall thickness can be 
seen in Figure 10. Coolant speed was held constant at 25 
m/ s . Coolant bulk temperature was held constant at 300 K. 
Recall from the analysis in Chapter 3 that boiling must not 
be permitted. For a wall thickness of 1.5 mm a coolant 
pressure greater than 2000 psi is needed. Furthermore, for 
reasons stated earlier, we would like to have the wall as 
thick as possible. For a wall thickness of 2.5 mm a coolant 
pressure of 1500 psi is sufficient, but notice that the air 
side wall temperature is close to the corrosion temperature. 
The wall thickness of 2.0 mm provides a Margin of Safety for 
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Temperature of 50 °C, and a Factor of Safety for Stress of 
about 1.4. A coolant pressure of 2000 psi is sufficient to 
prevent boiling. Based on the information in Figure 10 a 
nozzle liner thickness of 2.0 mm was selected. 

Figure 11 shows the effect of varying the coolant 
speed. In this case the liner wall thickness was held con- 
stant at 2.0 mm, and the coolant bulk temperature was held 
constant at 300 K. Again, the boiling temperatures of water 
at various pressures are superimposed. A coolant speed of 
20 m/s results in a coolant side wall temperature which can- 
not prevent boiling for water pressure at 2000 psia. A 
coolant speed of 30 m/ s causes a significant reduction in 
coolant side wall temperature, but does not significantly 
affect the Factor of Safety for Stress. Figure 11 indicates 
that a coolant speed of 25 m/s and a coolant pressure of 
2000 psia is optimum. 

Figure 12 shows the effect of varying the coolant bulk 
temperature. In this case we see that the bulk temperature 
over a range of 300 to 400 Kelvin has little effect on the 
wall temperatures or the Factor of Safety for Stress. 

Figure 13 shows the effect of using materials other 
than TZM. From this figure we see that TZM is the only 
material with Factor of Safety greater than 1.0 for liner 
thickness of 2 mm, coolant bulk temperature of 300 K, and 
coolant speed of 25 m/s. There is a possibility that some 
other set of conditions might result in greater Factors of 






15 


to 


coolant spee< 


rhe effect of co< 
L thickness 2*0 : 








Safety for the other materials. However, experience has 
shown that the molybdenum alloy TZM always outperforms the 
other candidate materials. 
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Start-up Considerations 

Next we consider start-up conditions. We define start- 
up conditions as wind tunnel operating conditions below the 
design 2000 psia pressure and 3460 °R temperature. The 
results of this analysis will predict the nozzle wall 
temperature over a range of start-up conditions. These 
predictions can be compared with temperature measurements 
made during start-up operation. If measured temperatures 
match predicted temperatures, we then can be confident that 
the nozzle will perform well at 3460 °R and 2000 psia. 

Analysis of the final nozzle design were performed at 
the start-up temperatures and pressures shown in Table 4. 

The result of this work appears in Figures 14, 15, 16 and 
17. In Figures 14, 15 and 16 the coolant wall temperatures 
are presented over the length of the nozzle for inlet air 
pressures of 1000, 1500, 2000 psia and inlet air tempera- 
tures of 2300 °R, 3200 °R and 3460 °R. Temperature 
measurements may be made with a thermocouple located away 
from the critical area near the throat. Such a location 
would be 2 cm downstream of the throat. Figure 17 predicts 
the coolant side wall temperature at 2.0 cm downstream of 
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the throat during steady state operation at inlet air 
pressures of 1000, 1500, and 2000 psia, and over a range of 
inlet air temperatures from 2300 °R to 3460 °R. 


Table 4. 


Start-up Temperatures and Pressures. 


Inlet Air Temperature 

Inlet Air Pressure 

2300 ° R 

1000 psia 

2300 

1500 

2300 

2000 

3200 

1000 

3200 

1500 

3200 

2000 

3460 

1000 

3460 

1500 
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Figure 15. Nozzle temperatures during start-up conditions. 











CHAPTER 5 


Conclusions and Recommendations 

The results of the analysis support the following con- 
clusions with regard to the nozzle shown in Figure 6: 

1. The material should be TZM as specified in the 
Metals Handbook [9] . This material can be obtained from 

| AMAX Specialty Metals Corp. [6]. Viking Metallurgical [7] 

or Northwest Industries [8] appear capable of fabricating 

I 

! the design. 

i 

I 

2. The wall thickness should be 2.0 mm along the 
length of the wall adjacent to the coolant passage surround- 
ing the throat area. The wall then tapers to 4.0 mm thick 
at the downstream flange, and 4.0 mm thick at the upstream 
end . 

3. The coolant should be water pressurized to 2000 
psi, with a bulk temperature around 300 K and flowing at 25 
m/s in the coolant passages. This corresponds to a flow 
rate of 127 gal/min. 

We recommend the installation of a thermocouple 2.0 cm 
downstream of the throat on the coolant side of the nozzle 
wall. Temperature measurements from the thermocouple may be 
compared with temperatures predicted from the analysis pro- 
viding a means of validating the analysis. 
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APPENDIX A 


Optimization of fin geometry 


The differential equation describing steady state heat 
transfer in a radial fin with rectangular profile is given 
by Kern and Kraus [9] as: 


r 2^e + r d9 

dr^ dr 


2 h 
k L 


r 2 0 = 0 


(Al) 


where 9 is at, and is the fin thickness as shown in 

Figure Al . We let m=(2h / k ) ^ , and note that Equation Al 
is Bessel's modified differential equation with general 
solution : 


0 = C ± I Q (mr) + C 2 K q ( mr) 


Substituting boundary conditions 


r = r 


r = r . 


e = e, 


d9 _ 


dr 


= 0 


into Equation A2 leads to 


S 0 = C 1 ^'“o 1 + C 2 K 0 ( "" r o> 
0 = I 1 (mr e ) + C 2 K 1 (mr e ) 


This leads to 


a 0 O [ K 1 (n»r e ) I Q (mr ) + I^mr^ K Q (mr) 
^ I 0 (mr o ) K i (mr e ) + (m O K n( m O 


0 


( A2) 


(A3a) 


(A3b) 


(A4a) 
( A4b) 


(A5) 
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Using Fourier's law of heat conduction at the base of 
the fin we have 

q = - 2U k r o ko • (A6) 

r=r o 

Differentiating Equation A5 at r=r Q and substituting into 
Equation A6 we have 

a. . = 21Tk r ^ a 

H with fin O °0 v ollQ(mr 0 ) K^(mr e ) + I^Cmrg) K Q (mc 0 ) 

(A7) 

We now calculate the heat transfer that would result if 
there were no fin: 


^without fin 2 ^ o r < 


Comparing the heat transfer obtained with and without the 
fin : 


-l 


^w/o 

f in 

- 2l>r o io 6. 

h 

(*r e ) 

K l (mr o ) “ K i< mr e ) I i (mr o ) 

%/ 

f in 

2 It it r <S _ m 

o v o 

*7 

_ l 0^ r o ) 

K^mrg) + I^mr^ K Q (mr 0 )_ 


2 

but m = (2h / ) . After cancellation: 
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The term in brackets consists of four Bessel function coef- 
ficients, which are transcendental functions. They depend 
on m and the size of the fin. For the magnitudes encoun- 
tered in the present study the term in brackets varies from 
0.7 to 0.9. We will assign the conservative value of 0.7 to 
this term. Clearly, we want the ratio of heat transfer 
shown on the left hand side of Equation A9 to be less than 
one. However, in order to arrive at a general guideline in 
the optimization process we will temporarily set Equation A9 
equal to unity. Then, carrying out the algebra and solving 
for the Biot number: 


Bi = (h 8 c )/k * 1.0 

Thus, in order to justify the use of fins, this ratio 
would need to be significantly less than 1.0. Note that as 
the Biot number decreases, cooling fins become more and more 
justifiable. 

Substituting typical values from our study for the 
variables we have 

0 O = 222 Kelvin 

^> o = 0 .002 meters 

r Q = 0.00769 meters 

r = 0.0102 meters 
e 

h = 100,000 W/m 2 . K 


k 


100 W/m-K 
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From Jahnke, Emde and Losch [10] we have 


1 0 {mt o ) 

= 

I 0 (7.7) 

=r 

323.1 

I l (mr o ) 

= 

I-L ( 7 . 7) 

= 

301.1 

*1 ( mr e ) 

= 

Ij^dO) 

s 

2671 

K 0 (mr o ) 

= 

K q (7.7) 

= 

20142 

K l (mr o ] 

= 

K x ( 7 .7) 

= 

21412 

K 1 ( mr e ) 

= 

K x (10) 

= 

18648 


Substituting these values into Equation A9 we have: 

(h J 0 )/k = 1.16 . (A10) 

This ratio indicates that the presence of a fin actu- 
ally decreases the amount of heat transfer. Due to the high 
convection coefficient on the coolant side as compared to 
the thermal conductivity and thickness of the fins we con- 
clude that in this case cooling fins cannot improve heat 


transfer . 



APPENDIX B 


Air Side Heat Transfer 

In a hypersonic wind tunnel acceleration of the gas (in 
this case dry air) is so great that simplifying assumptions 
usually made for duct flow no longer apply. Bartz [1] 
points out that the axial pressure gradient terms in the 
momentum and energy equations must be accounted for. 

Using the integral momentum and energy equations we can 
solve for the Stanton number, St, as a function of axial 
location. The air side heat transfer coefficient as a func- 
tion of axial location will then be given by 

h a = (St C P U e ) alr , (Bl) 

where U is the local free stream speed of the air, and 
and Cp are properties of air evaluated at the local free 
stream conditions. Bartz developed the equations necessary 
to solve for the local air side heat transfer coefficient. 
The assumptions made for this development are reproduced 
here in their entirety. 

1. The flow is axisymmetric and steady without tan- 
gential components of velocity. 

2. The boundary layer is confined to a distance from 
the wall which is small compared with the distance from 
the axis of symmetry. 
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3. The only forces acting on the gas are those due to 
pressure gradients and to skin friction at the wall. 

4. The only changes in total enthalpy in the flow dir- 
ection are those due to heat flux through the wall. 

5. The flow immediately outside the boundary layer is 
reversible and adiabatic and parallel to the wall. 

6. Static pressure is constant through the boundary 
layer perpendicular to the wall. 

7. The gas is perfect; however, the restriction that 
specific heats be constant can be removed in computing 
the driving potential for heat flux. 

8. The gas has a constant Prandtl number, a viscosity 
which varies as a power of the temperature, and a con- 
stant adiabatic recovery factor. 

9. The skin-friction coefficient is the same as for 
constant-pressure constant-wall-temperature flow on a 
flat plate at the same free-stream conditions, wall 
temperature, and momentum thickness. 

10. The Stanton number is the same as for constant- 
pressure constant-wall-temperature flow on a flat plate 
at the same free-stream conditions, wall temperature, 
and momentum thickness. 

11. The Stanton number for unequal momentum and energy 
thicknesses is that for equal thicknesses multiplied by 
( M/ 0 ) n , where n is a small "interaction exponent". 
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12. Heat transfer affects the skin-friction coeffi- 
cient in either one of two ways: (a) there is no 

effect, and is the same as for adiabatic flow, or 

(b) is the same as for adiabatic incompressible flow 
at a density and viscosity evaluated at the arithmetic 
mean between the actual wall temperature and the free- 
stream static temperature. 

13. The Stanton number for equal momentum and energy 
thicknesses is related to the skin-friction coefficient 
by von Harman' s form of Reynolds' analogy. 

14. Any chemical reactions in the boundary layer 
affect only the driving potential for heat flux. 

15. The boundary-layer shape parameters 8/ <S , a/£ , 

and &*/ 9 are those of 1/7 power profiles of velocity 
and of the difference between stagnation and wall tem- 
perature. Such profiles are typical of turbulent 
boundary layers on flat plates. 

16. Heat transfer by thermal radiation is negligible 
compared with convection. 

17. There is no significant net mass transfer from 
wall to gas or gas to wall. 


The applicable equations that result from Bartz's work 
involve the axial gradients of the momentum and energy 


thicknesses : 
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With initial values for 9 and we can set up a 

finite difference numerical scheme for the solution of Equa- 
tions B1 through B5 as a’ function of axial location. Bartz 
showed that values for 0 and at the throat are quite 

independent of the initial values. 

It should be noted that both Equations B2 and B3 con- 
tain the axial gradient of Mach number. This gradient is 
relatively small in the upstream section and the axial step 
size is allowed to be correspondingly large (0.25 cm in this 
case). However, just downstream of the throat this gradient 
becomes quite large and a much smaller step size is required 
(0.125 mm) . The smaller step size means more time is 
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required for solution. Perhaps more importantly the possi- 
bility for cumulative round-off error increases. 

Two other methods for calculation of the air side heat 
transfer coefficient were considered. The first was a sim- 
plified equation from Bartz [12] meant to be used for quick 
estimates of h . However, the results were shown to be 

cl 

always higher than Bartz's more involved solution (see ref- 
erence [1]). Due to the critical nature of conditions at 
the throat, it was felt that an artificially high h would 
impose unnecessary restrictions on the design. For example, 
an extremely high h_ might force the design to have a very 
thin liner wall thickness at the throat, perhaps around one 
millimeter. Such a thin wall might compromise the design in 
terms of handling during installation, or erosive "sand- 
blasting" during operation. The simple calculation for h 

3 

was rejected in favor of a more accurate model. 

The second alternative method for calculation of the 
air side heat transfer coefficient was given in Rohsenow and 
Hartnett [2] . Here, a number of empirical formulas were 
used. This method yielded an h within about 30% of the 
other methods, although it was consistently low. Again, the 
accuracy of the results might be questioned in terms of the 
formula's applicability to hypersonic flow. 


APPENDIX C 


Coolant Side Heat Transfer 

The cooling system configuration shown in Figure 6 of 

i 

the main text consists of three coolant passages. Each pas- 
sage is supplied by coolant transitioning from a circular 
tube to a rectangular passage. From the rectangular passage 
the flow is split and diverted over and around the nozzle 

f 

I 

liner as two curved rectangular ducts. This is shown in 
Figure Cl. 

Our objective is to predict, using published litera- 
ture, the heat transfer from the nozzle liner to the coolant 
in the curved rectangular passageways. A problem exists, 
however, in that the literature does not directly give heat 
transfer information for our particular application. 

Instead, the heat transfer must deduced based on information 
available for other similar flow passage configurations and 
flow conditions. 

Our cooling system falls under the classification of 
internal duct flow. As a start, for flows that are char- 
acterized by large property variations, the following 
equation is recommended by Sieder and Tate [13] 

Nu d = 0.027 Re D 4/5 Pr 1/3 (/*//!„) °* 14 f (Cl) 
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with the following restrictions: 

circular duct flow 
turbulent flow 
fully developed flow 
smooth wall 

uniform circumferential heat flux 

0. 7. < Pr < 16,700 

Re D > 10,000 

L/D > 10 

fluid properties at coolant bulk temperature 
(except yU w ) 

Our flow passage and flow conditions do not conform to 
the restrictions imposed by Equation Cl. The deviations are 

1. rectangular cross section duct 

2. flow is not fully developed 

3. circumferential heat flux variation 

4. curved flow passage 

Each one of the above deviations are now discussed with 
the corresponding effect on Equation Cl. 

1. Rectangular cross section duct. Incropera and De Witt 
[14] report that Equation Cl may be used for noncircular 
ducts if the hydraulic diameter is used in the Nusselt and 
Reynolds numbers. The hydraulic diameter is defined as 
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D h = 4A C /P • (C2) 

This substitution is confirmed by Gambill and Bundy 

[15] as reported in Rohsenow and Hartnett [2], where Equa- 
tion Cl fit the data very well for water in a rounded corner 
rectangular channel. 

2. Flow is not fully developed. The flow is definitely 
developing thermally (start of the thermal boundary layer) 
at the point where it splits to travel around the nozzle 
liner. Also, at this point there is some argument that it 
may be hydrodynamically developing (start of the velocity 
boundary layer) as well. For the combined thermal and 
hydrodynamic entry length situation considerable data are 
available for circular tubes [16,17,18,19,20], all of which 
can be summarized by 

^ U D entrance ^ ^ U D fully developed . (d) 

region region 

Since our configuration is not tube flow and clearly not a 
well defined entrance region we will not attempt to account 
for this effect. Any error associated with this omission is 
on the conservative side. That is, the heat transfer coef- 
ficient will be lower as a result of not accounting for de- 
veloping flow. 

3. Circumferential heat flux variations. The coolant pas- 
sage is bounded on the top by the structural housing and on 
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the bottom by the nozzle liner as shown in Figure C2 . The 
nozzle liner is heated, while the structural housing is not. 
This results in what is called a "circumferential heat flux 
variation". Rohsenow [2] reports the work of Reynolds [21], 
and Sparrow and Lin [22] that applies to circular tubes with 
fully developed flow under circumferential heat flux varia- 
tions. We will not attempt to account for this effect due 
to differences in our flow and flow channel geometry from 
those reported. Again, any error associated with this omis- 
sion is on the conservative side. The thermal boundary 
layer is building up on only one side of the channel. This 
has the effect of a longer thermal entrance region with 
associated higher heat transfer rates. 

4. Curved flow passage. A number of investigations [23, 
24,25] examined heat transfer and pressure drop in curved 
passages. From these papers a general conclusion can be 
drawn that heat transfer is greater on the concave side of 
the flow passage compared to the convex side. This is un- 
fortunate since we want to cool the convex side only in our 
flow passage. Kreith [24] reports the results of his ex- 

4 

perimental work with water at Re Q = 5 X 10 m a curved 
parallel plate channel with r^/r 0 = This Reynolds 

number is close to our design Reynolds number. Addition- 
ally, our design radius ratio is r^/r = 0.7. Kreith gives 

Nu /Nu as a function of Prandtl number. Our T, 

concave 7 convex b 
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will range 
number and 


from 300 K to 350 K, 

Nu /Nu 

concave 7 convex 


which corresponds to Prandtl 
as shown in Table Cl. 


Table Cl. 

Heat Transfer in a Curved Parallel Plate Passage. 



Pr 

Nu 

concave 

T b 

Nu 

convex 

300 

5.85 

1.55 

350 

2.33 

1.49 


In Table 1 we will take an average Nusselt number ratio 
of 1.5. From Yang and Liao [25], who performed experiments 
in rectangular ducts with 180 bend, we obtain the general 
relationship 


(Nu 


concave ^ u straight^ ^ u straight 


- Nu 


convex 


(C4) 


Recallinq that Nu /Nu = 1.5 gives 

y u concave 7 convex ^ 


^convex ~ ^straight' 


(C5) 


Nu . . .. is the same as the Nu^ in Equation 1. 

straight d 

Therefore Equation 1 is substituted into Equation 5 to give 
the final result for heat transfer on the coolant side 


h c = k/D h °* 0216 Re ° ’ 8 Pr 1/3 ( yU / yU w ) 0,14 


(C6) 
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Figure Cl. Coolant passage at nozzle throat. 



APPENDIX D 


Numerical Methods of Solution 
Flow Analysis. Local Mach number is an implicit func- 
tion of local area divided by throat area: 


_A_ _ 1 (1 + 0.2 Ma 2 ) 3 

A* * Ma 1.728 


(Dl) 


Equation Dl applies to air with # = 1.4. Local area is a 

2 

function of local radius using A = ft r . Local radius a 
function of axial distance is available to the FORTRAN code 
as a subroutine. Mach number as a function of axial distance 
is therefore obtained by solving Equation Dl using a numeri- 
cal regula falsi technique. With a solution for local Mach 
number, local free stream temperature and pressure are avail- 
able from compressible flow theory: 


T = 
■ L oo 


1 + 0.2 Ma' 


(D2) 


Poo = P Q (1 + 0.2 Ma 2 ) -3 * 5 


(D3) 


Thermal Analysis. A method of successive substitution 
is required for the thermal analysis in the cooled portion of 

the nozzle. An initial value is obtained for T w , c usin 9 
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Equation 8 of the main text and initial guesses for h c and 
T . The guess for h can now be updated using Equation 6 

W g 3, C 

from the main text. Furthermore, T can now be solved for 

w , a 

using the new estimate for T and a guessed value for h in 

Equation 7. If the new estimate for T a is significantly 

different from the current estimate, then it becomes the 

current estimate. The process is repeated starting with a 

solution for T using the new estimate for h . The new 

estimate for T is significantly different form the 

w , a 

current estimate if it differs by more than one degree 
Celsius . 

Note that this successive substitution scheme does not 
modify h . Once T has converged we begin another loop. 

The air side wall temperature is now used to calculate a new 
value of h with Equation 5. If this new h_ differs from the 

3 3 

2 

current value by more than one watt/m .°C, then it becomes 
the current value. The loop that solves for T a is repeated 
first; then h is retested. When h finally converges we 

3 a 

have the air side wall temperature and heat transfer coef- 
ficient, and the coolant side wall temperature and heat 
transfer coefficient as a function of axial position. 

Stress Analysis. Equations 10 through 19 in the main 
text deal with the stress analysis. The AT used in Equation 
10 and Equation 12 is the difference between air and coolant 
side wall temperatures. We have those temperatures from the 
thermal analysis. Also needed is the local free stream air 
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pressure. This was obtained fron. the flow analysis. The 
stress analysis is therefore a relatively straight forward 
application of Equations 10 through 19. 

The FORTRAN code flow chart is given in Figure Dl. The 
code listing is given at the end of this appendix. 



Figure D1 . 


FORTRAN Code Flow Chart 
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FORTRAN CODE 


PROGRAM MOD2 


CHARACTER* 100 I DENT 

CHARACTER * 50 AGA I N , HALT , RES , METAL , UNITS . TYPE , CYL 
EXTERNAL FMA , TOP, BOT 

COMMON/FT/ST , TR , TW , TO , DRDZ , MA , DELMA , STEP , RAD , CF , DSDTH, T I NF , DELTW 
COMMON/AR/AR 

COMMON/COEF/DELP. DH. L, RO, NU, FRIC. R1 , XI , X3 

COMMON/DF/DF 

INTEGER STATN, POS, NEG 

REAL MDOT , NU,MA,L,K,MUO,MUAIR,NUM, NTHTA . NPHI 

REAL TWAC-100: 85) , TWCC - 100 : 85) , YLDl C -100 : 85}, YLD20100: 85) 

REAL HAC' 100 : 85) , HCC ” 100 : 85) , QC~ 100 : 853 , QPPC-100: 85) 

REAL EACO: 15), AACO: 15), PACO: 15), AKCO: 15) , YACO: 15) 

REAL THC10),UMNC10),TBKC10) 

REAL Kl,K2,K3,K4 
0PENC6,FILE«'PRN: ' ) 

DF-400.0 

* " THIS DIVIDING FACTOR IS USED TO LESSEN THE ROUND OFF ERROR 

* IN FITTING A POLYNOMIAL CURVE TO THE FUNCTIONAL RELATION 

* BETWEEN THE MATERIAL PROPERTIES AND TEMPERATURE 

DO 3 1-0, 15 
EACD-0.0 
AACD-0.0 
PA( I )»0 • 0 
AKCD-0.0 
YAC I )=0 . 0 
3 CONTINUE 

555 CONTINUE 

CALL DATA INC EA , AA , PA , AK , YA , METAL) 

PC-13. 7895E6 

* ~ PA — COOLANT PRESSURE C2000 PSI) 

554 CONTINUE 

4321 CONTINUE 

PRINT* , ' IDENTIFIER: * 

READ* C A) ' , IDENT 

DISTC-1.0 
TMPC-i.O 
VELC-1.0 
HCNV-1 . 0 

DNSTR-0 . 12 
DNTHK-0 . 004 
UPSTR-0 . 10 
UPTHK-0 . 004 
X3-0.004 
X2-0.08 

T0-3460. 0 
TO-TO/1 .8 
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P0-2000.0 
PO=PO * 6894 . 8 

TCRIT-760.0 

TCRIT-TCRIT+273.0 

5432 PRINT*.' ENTER M S*‘ FOR SINGLE RUN OR ,# P # ’ FOR PARAMETRIC OUTPUT ' 
READ' C A) * , TYPE 
IF CTYPE . EQ. 'P' DTHEN 

PRINT*,' NUMBER OF THROAT THICKNESSES TO BE TESTED: ' 

READ * , N1 

IF CUNITS .NE. 'EN' )THEN 

PRINT*,' ENTER THICKNESS IN METERS ' 

ELSE 

PRINT*,' ENTER THICKNESS IN INCHES ' 

END IF 

DO 31 I-l.Nl 

PRINT*,' ENTER THICKNESS' , I, ' : ' 

31 READ* , THC I 3 

PRINT*,' NUMBER OF COOLANT SPEEDS TO BE TESTED: 

READ * , N2 

IF CUNITS .NE. ’ EN' )THEN 

PRINT*, ' ENTER SPEED IN METERS/SECOND ' 

ELSE 

PRINT*,' ENTER SPEED IN FEET/SECOND ' 

END IF 

DO 32 I-1.N2 

PRINT* , ' ENTER SPEED' , I , ' : ' 

32 READ* , UMNC I !> 

PRINT*,' NUMBER OF COOLANT BULK TEMPERATURES TO BE TESTED: 

READ * , N3 

IF CUNITS .NE. 'EN'DTHEN 

PRINT*,' ENTER TEMPERATURE IN KELVIN ' 

ELSE 

PRINT*,' ENTER TEMPERATURE IN RANKINE ' 

END IF 

DO 33 I-1.N3 

PRINT* , ' ENTER TEMP. ' , I , ' : ' 

33 READ * , TBKC I ) 

PRINT*,' REPEAT INPUT CYCLE CY/N3? ' 

READ' C AD ' , AGAIN 
I FC AGAIN .EQ. 'Y')GO TO 554 

ELSE 

PRINT*,' LINER THICKNESS AT THROAT' 

IF CUNITS .NE. 'EN'DTHEN 

PRINT*,' ENTER THICKNESS IN METERS: ' 

ELSE 

PRINT*,' ENTER THICKNESS IN INCHES: ' 

END IF 

READ* ,THTHK 
THTHK-THTHK* DISTC 


ORIGINAL PAGE IS 
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PRINT*.’ AVERAGE COOLANT SPEED AT THE THROAT’ 
IF CUNITS . NE. ’ EN’ 5THEN 

PRINT*.’ ENTER SPEED IN METERS/SECOND: ’ 

ELSE 

PRINT*.' ENTER SPEED IN FEET/SECOND: ' 

END IF 
READ* ,UMT 
UMT-UMT*VELC 

PRINT*.’ COOLANT BULK TEMPERATURE’ 

IF CUNITS .NE. 'EN’OTHEN 

PRINT*.’ ENTER TEMPERATURE IN KELVIN: ’ 

ELSE 

PRINT* , ' ENTER TEMPERATURE IN RANKINE: ' 

END IF 
READ* , TB 
TB-TB* TMPC 

PRINT* . ' REPEAT INPUT CYCLE CY/ND? ' 

READ' CAD ’ , AGAIN 

IFCAGAIN . EQ. 'Y'DGO TO 554 

GO TO 6543 

ENDIF 

DO 20 IT-1, N1 
THTHK-THCITD 
DO 21 IB-1, N3 
TB-TBKC IBD 
DO 22 I U-l , N2 
UMT-UMNC I UD 

6543 CONTINUE 

DO 5 1 — 100,85 
TWACID-0.0 
TWCCID-0.0 
■YLDICO-O.O 
YLD2CID-0.0 
HACID-0.0 
HCCID-0.0 
5 CONTINUE 

TOL-O.O 

PI-3.1415926 

GAM-1.4 

RTH-O. 224* .0254 

* ^ METERS — THROAT DIA. 


STEP-0.0025 
RAD-RAD I USC-0 . 235D 
HAC-48D-1200.0 
HCC -48 D *50000 . 0 
TWAC-48 D-500 . 0 
TM-430.0 

AR-RAD* RAD/RTH/RTH 

Rl-RAD/COSC AT ANCC RAD-RAD I USC-0 . 105} D/ STEP D D 


CALL SOLVEC 1 • OE-9 . l.O.O.DOl, 500 , MA , FMA) 

T I NF-TO/C 1 .0+0.2*MA*MA) 

TR-TI NF+Q . 87* CTO-TINF) 

P I NF-PO*CT INF/TO)* *3.5 
UA I R-MA * 20 . 0 * T I NF **0.5 
TWAC-94)-TR 

EPS1-1.0 
EPS2-1.0 
EPS3-. 00001 
EPS4-1 . OE- 10 
EPS5-1.0E-5 

TBR-1 . 8 * TB 

PROPERTIES OF WATER OVER APPLICABLE RANGE 
CK-0 . 4483+0 . 05022 * LOGCTB-273 . 0) 

CRO- 1001 .5-0. 2217* CTB-273.0) 

CPR«133.4*CTB-273.0)**C-0. 9354} 

CMU-O. 0254 *CTB-273 .0)* *00-9828) 

SPC-4200 . 0 
NU-CMU/CRO 

L-P I * CRTH+THTHK+X3/2 . 0) 

AREA-X2 * X3 
PER-2 . 0 * CX2+X3) 

DH-4.0* AREA/PER 
RE-UMT* DH/NU 

MDOT-2 . 0*CR0* AREA*UMT 

THIS IS THE TOTAL MDOT CTOP AND BOTTOM CHANNELS) 
UPBGN-0.0254 
DNBGN-O. 0381 
PHI-0.00001 
THTA-O. 00001 

****** DO 10 INCREMENTS THE AXIAL POSITION IN THE NOZZLE 
DO 10 STATN— 93,4 

AX-STEP* STATN 

IF C STATN . LT . -16)THEN 

XI- . 004-C48+STATND/400. 0* C . Q04-THTHK)/0 . 08 
ELSE IF CSTATN . GT . 16)THEN 

Xl-THTHK+C STATN- 16)/ 400 . 0* CO . 004-THTHK)/0 . 08 
ELSE 

Xl-THTHK 
END IF 

ROLD-RAD 
RAD- RAD I USC AX) 

AR-RAD* RAD/RTH/RTH 
R1-RAD/C0SCATANCCRAD-R0LD)/STEP)) 

RM-R1+X1/2 . 0 
R2-R1+X1 
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DELR-RAD-ROLD 

DRDZ-DELR/STEP 

DELX-C STEP * *2+DELR* *2 VO. 5 

X-X+DELX 

OLDMA-MA 

IF CSTATN .GE. 0) THEN 

XA-1.0 

XC-16.0 

ELSE 

XA-0. 00001 

XC-1.0 

ENDIF 

CALL SOLVE CXA, XC, . 001 , 500, MA , FMA) 

TINF*T0/C1 . 0+0 . 2 * MA *MA) 

P I NF*PO*CT INF/TOD* *3.5 
TR-TINF+O . 87* CTO-TI NF) 

^ RECOVERY TEMP. — ALL UNITS ARE SI COMPATIBLE 
UAIR-MA*20.0*TINF* *0.5 
DELMA-MA-OLDMA 

ROA I R-P I NF/TI NF/287 . 0 

CPA I R-CO . 93608+0 . 00019834 * T 1 NF) * 1000 . 0 
MUA I R-CO . 05131 " T I NF * * 0 . 63684) *0 . 00001 
PRAIR-O. 78145-0. 01 169 *LOGCTINF) 

*** THERMAL ANALYSIS ******** 

I.FCSTATN . LT . -48)THEN 

TWACSTATND-TR 

TW-TR 

DELTW«TWACSTATN)-TWACSTATN-i) 

NUM-TRAPCTOP.O.O, 1.0, 100) 

DENOM-ROMBCBOT , 0 . 0 , 1 . 0 , 5 , EPS5 , 0) 

DSDTH-C 1.-7. *TINF*NUM/TWACSTATN))/C7. *TINF* DENOM/TWACSTATN)) 
RET “ROA I R* UAI R *THTA/MUA I R 

CF-0.0256*RET* *00. 25)“CC.5*CTWA(STATN)/TINF+1. 0))* *00.6) 

Ki-FTHTACTHTA) 

K2-FTHTACTHTA+0 . 5*STEP*K1) 

K3-FTHTACTHTA+Q . 5*STEP*K2) 

K4-FTHTACTHTA + STEP * K3) 

THTA-THTA+STEP/6 .0*CK1+2.G*K2+2.G*K3+K4) 

K1-F1PHI CPHI ) 

K2-F1PHI C PHI +0.5* STEP *K1) 

K3-F1PHICPHI+Q. 5*STEP* K2) 

K4-F1PHICPHI+STEP*K3) 

PHI -PHI + STEP/ 6.0*CK1+2.0*K2+2.0*K3+K4) 


ELSE 
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K-YYCTM, AKO 
15 CONTINUE 

TWCCSTATN0-CK*TWACSTATN0/LOGCR2/R10+HCCSTATN0*R2*TB0/ 

& CHCCSTATN0*R2+K/LOGCR2/R100 
TBR-1.8*TB 

YMU-0 . 0254 ■ CTWCCSTATNO -273 . 05 * * C-0 . 9828 0 

HCCSTATNO-O. 0216 *CK* CUM *CROO* *0. 8*CPR* *C . 333333330/CCMU* *0.8* 
& DH* *0. 25 *CCMU/YMU0* *0. 14 
RE-UM* DH/NU 

TM-CTWCCSTATNO+TWACSTATNOO/2 . 0 
K-YYCTM, AKO 

TWAN-CHACSTATN0*R1 * TR+K* TWCCSTATNO/ L0GCR2/R1 0 0/ 

& CK/LOGCR2/R10+HACSTATN0*R10 
IFCABSCTWAN-TWACSTATNOO .GT. EPS10THEN 
TWACSTATNO-TWAN 
GO TO 15 
END IF 

TWACSTATNO-TWAN 

DELTW-TWACSTATNO-TWACSTATN-IO 
TW-TWACSTATNO 
NUM-TRAPCTOP .0.0.1.0,1000 
DENOM-ROMBCBOT, 0. 0, 1.0,5, EPS5 , 0} 

DSDTH-C1.-7. *TINF*NUM/TWACSTATNOO/C7. *TINF*DENOM/TWACSTATNOO 
RET-ROAIR*UAIR*THTA/MUAIR 

CF-0.0256*RET**C-0.250*C0.5*CTWACSTATNO/TINF+1 .000** C-0. 60 

Kl-FTHTACTHTAO 
K2-FTHTACTHTA+0. 5*STEP*K10 
K3-FTHTACTHTA+0. 5*STEP*K20 
K4-FTHTACTHTA+STEP*K30 

THTA-THTA+STEP/6 . 0 * C Kl+2 . 0*K2+2 .O*K3+K40 

ST-CPHI/THTAO* *C-0. 1 0 *CF/2 . 0/ C 1 • 0-5 . 0-* CCF/2 . OO * *0 . 5* 

& Cl .0-PRAIR+L0GC6. 0/C5 . 0* PRAI R+l .00000 

K1-F2PHI CPHI 0 
K2-F2PHI CPHI +0.5 'STEP* K10 
K3-F2PHI CPHI +0.5* STEP *K20 
K4-F2PHI CPHI +STEP*K3 0 

PHI -PHI +STEP/6 . 0* CK1+2 . 0* K2+2 . O*K3 + K40 

HAN-ST* ROAIR'CPAIR* UA IR 

IFCABSCHAN-HACSTATNOO ■ GT . EPS20 THEN 
HACSTATNO-HAN 

GO TO 15 
END IF 
END IF 


HACSTATNO-HAN 
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deltat-twacstatnd-twccstatnd 
TWCCSTATN+ 1 5-TWCCSTATND 
TWACSTATN+n-TWACSTATN) 

HCCSTATN+n-HCCSTATND 
HACSTATN+1 D-HACSTATN} 

•**•* STRESS ANALYSIS * 

DR-R2 * R2-R1 * R1 

F1R2-YYCTWACSTATND . AAD * YYCTWACSTATN5 , EA) ‘DELTAT* 
& C1.-2. *R2*R2*LOGCR2/Rn/DRV 
& C2 . * C 1 • -YYCTWACSTATND .PAD} * LOGCR2/R1 

STH1-F1R2+PINF*CR2*R2+R1*R1 D/DR-2 . *PC*R2‘R2/DR 

SZ1-F1R2 

SRI— PINF 

FlRl-YYCTWACSTATND, AAD * YYCTWACSTATND . EA 5 * DELTAT * 
& Cl - -2. * R1 * R1 * LOGCR2/R1 3/DRD/ 

& C2 . *C1.-YYCTWACSTATN3,PAD3*L0GCR2/R1)) 

STH2— F1R1+2. 0*PINF*R1*R1/DR-PC*CR2*R2+R1*R1)/DR 
SZ2— F1R1 
SR2— PC 


YLD1CSTATND-C . 5 * C CSTH1-SR1 5 * ‘2+CSTH1-SZ1 D * ‘2+CSZ1-SR1 } * *2DD * * . 5 
YLD2CSTATND-C . 5 * CCSTH2-SR2D * ‘2+CSTH2-SZ2D* * 2 + CSZ2-SR2D * *2}} * * .5 

LOCAL HEAT TRANSFER Q FOR COOLED PORTION OF NOZZLE 

QCSTATN}-2.0*P! *RAD‘ ABSCSTEPD * HCCSTATND * CTWCCSTATND-TB} 

qppcstatnd-hccstatnvctwccstatnd-tbd 

........ L00 p 10 


10 CONTINUE 

STEP-0.000125 

........... D0 n 

DO 11 STATN-5 , 85 

AX-STEP*STATN+0. 01 

Xl-THTHK 

UM-UMT 

ROLD-RAD 
RAD-RAD I USC AX D 
AR-RAD* RAD/RTH/RTH 

Rl-RAD/COSC AT ANC CRAD-ROLD D/STEP) } 

RM-R1+X1/2 - 0 

R2-R1+X1 


DELR-RAD-ROLD 
DRDZ-DELR/STEP 
DELX-CSTEP* *2 + DELR* *20* *0. 5 
X-X+DELX 

OLDMA-MA 

IF CSTATN .GE. 03 THEN 

XA-1.0 

XC-16.0 

ELSE 

XA-O. 00001 
XOl.O 
END IF 

CALL SOLVE CXA , XC . . 001 , 500 , MA , FMA3 

TINF-T0/C1.0+0.2*MA*MA3 
PINF-PO* CTI NF/T03 **3.5 
TR-TINF+0.87*CT0-TINF3 

~ RECOVERY TEMP. -- ALL UNITS ARE SI COMPATIBLE 
UAIR=MA*20. O’TINF* *0.5 
DELMA-MA-OLDMA 

ROA I R-P I NF/T I NF/287 . 0 

CPA I R-CO . 93608+0 . 00019834 * T INFO * 1000 . 0 
MU A IR-CO. 05131 *T INF* *0.636840*0.00001 
PRA I R=0. 78145-0. 01 169* LOGCT INFO 

...... THERMAL ANALYSIS (110 ******** 

K-YYCTM, AKO 
16 CONTINUE 

TWCCSTATNO-CK*TWACSTATNO/LOGCR2/R10+HCCSTATNO’R2*TBO/ 

6 CHCCSTATNO*R2+K/LOGCR2/R130 
TBR-1 . 8*TB 

YMU-O. 0254 ‘CTWCCSTATNO-273 .03’ *C-0. 9828 0 

HCCSTATNO-O. 0216* CK* CUM* CR03* *0.8*CPR* *C • 333333330/CCMU * *0 . 8 * 
& DH* *0. 2 3 * CCMU/ YMU3 * *0.14 
RE-UM * DH/NU 

TM-CTWCCSTATNO+TWACSTATNOO/2 . 0 
K-YYCTM, AKO 

TWAN-CHACSTATN0*Rl*TR+K*TWCCSTATN0/LOGCR2/R100/ 

& CK/L0GCR2/R1 0+HACSTATN3 *R1 0 
IFCABSCTWAN-TWACSTATNOO .GT. EPS10THEN 
TWACSTATNO-TWAN 
GO TO 16 
END IF 

TWACSTATNO-TWAN 
DELTW-TWACSTATN3-TWACSTATN-10 
TW-TWACSTATNO 
NUM-TRAPCTOP, 0.0,1 .0,1000 
DENOM-ROMBCBOT ,0.0, 1.0,5, EPS5 , OO 

DSDTH-C1 • -7. *TINF’NUM/TWACSTATN00/C7. *TINF*DENOM/TWACSTATN33 
RET-ROAIR*UAIR*THTA/MUAIR 

CF-O. 02 56* RET* *C~O.250*CO. 5* CTWACSTATNO/TINF+1 . OOO * ’ C-0 . 60 
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ki-fthtacthta:> 

K2-FTHTACTHTA+0 . 5 4 STEP 4 Ki ) 

K3 -FTHT AC THT A+O . 5 4 STEP 4 K2 3 
K4-FTHTACTHTA+STEP-K3) 

THTA-THTA+STEP/6.0-CK1+2. O 4 K2+2 . 0 4 K3+K4) 

ST-CPHI/THTAD 4 4 C-0. 1 ) 4 CF/2 . 0/ C 1 . 0-5 . 0 4 CCF/2 . 0) 4 ‘0 . 5* 

& C1.0-PRAIR+LOGC6.0/C5.0 4 PRAIR+1.0)3)D 

K1-F2PHICPHI) 

K2-F2PHICPHI+0. 5 4 STEP 4 KID 
K3-F2PHI CPHI +0 . 5 4 STEP 4 K2) 

K4-F2PHI CPHI+STEP 4 K3} 

PHI ~PHI +STEP/6 . 0 4 CK1+2 . 0* K2+2 . 0 4 K3 + K4) 

HAN-ST 4 ROAIR 4 CPAIR 4 UAIR 

IFCABSCHAN-HACSTATN)) . GT . EPS2) THEN 
HACSTATNO-HAN 

GO TO 16 
ENDIF 

HACSTATNO-HAN 

DELTAT*TWACSTATN])-TWCCSTATN) 

TWCCSTATN+1>TWCCSTATN) 

TWACSTATN+D-TWACSTATND 

HCCSTATN+D-HCCSTATND 

HACSTATN+1D-HACSTATND 

44444 STRESS ANALYSIS ********** 

DR=R2 * R2-R1 4 R1 

F1R2~Y YCTWACSTATN) , AAO 4 YYCTWACSTATN} , EA) 4 BELT AT 4 
& Cl- -2. 4 R2 4 R2 4 L0GCR2/R1)/DR)/ 

& C2. 4 C1.-YYCTWACSTATN),PAD: 4 LOGCR2/R10) 

STH1*F1R2+P I NF 4 CR2 4 R2+R1 4 R1)/DR~2. 4 PC 4 R2 4 R2/DR 

SZ1-F1R2 

SRI — PINF 

F1R1-YYCTWACSTATN} , AA) 4 YYCTWACSTATN3 , EA} 4 DELTAT 4 
6 C1.-2. 4 R1 4 R1 4 L0GCR2/R13/DR^/ 

& C2 . 4 C1.-YYCTWACSTATN).PAD^ 4 L0GCR2/R1D) 

STH2— FlRl+2.0 4 PINF 4 Ri 4 Rl/DR-PC 4 CR2 4 R2 + Ri 4 RlI)/DR 
SZ2 — F1R1 
SR2— PC. 

YLD1CSTATND-C.5 4 CCSTH1-SR1D 4 4 2+CSTH1-SZ 1 ) 4 ‘2 + CS21-SR1 ) 4 4 2) 3 4 4 . 5 
YLD2CSTATN3-C . 5 4 C CSTH2-SR2 3 4 4 2+CSTH2-SZ2 3 4 4 2 + C SZ2-SR2 3 4 4 2 0 3 4 4 . 5 

4 LOCAL HEAT TRANSFER Q FOR COOLED PORTION OF NOZZLE 
QCSTATN>2.0 4 PI 4 RAD 4 ABSC STEP} 4 HCCSTATN3 4 CTWCCSTATN3-TB3 
QPPCSTATN3~HCCST ATN3 * CTWC(STATN3~TB3 

********* LOOP 11 ************** 


11 CONTINUE 
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......... EQUATIONS MAKING USE OF THE LOCAL HEAT TRANSFER Q 

WRITEC6.80COIDENT 

800 F0RMATC5X, A/D 

WRI TEC6 , 801 DMETAL 

801 F0RMATC5X, ' MATERIAL IS ',A//D 
............ BLOCK J *■**•*■*•***• 

WRI TEC6 * 5020 

502 FORMAT C' LOG TWA TWC FSA FSC HA',7X.’HC' 

& * 8X i ' QPP MSA' D 

DO 57 1—28,4 
XL0C-I/4.0 
TMAR-TCRIT-TWACID 
FN1-YYCTWACID. YAD 
FD1-YLD1CO 
FF1-FN1/FD1 
FN2-YYCTWCCI0. YAO 
FD2-YLD2CID 
FF2-FN2/FD2 

IFCI .GE. -5 .AND. I .LE. 4DTHEN 

PRINT* , ' SY* ' , FN1 

PRINT*, * LOCAL STRESS- ',FD1 

PRINT 4 , ' F . S . - ' , FN1/FD1 

END IF 

WRITEC6, 500DXL0C, TWAC I D/TMPC , TWCC I D/TMPC , FF1 , 

& FF2 , HAC I 0 * HCNV , HCC I D * HCNV , QPPC I D . TM AR 
14 • ENGLISH UNITS FOR HTXFR COEFF ARE BTU/ CHR * F 4 FT2 0 

57 CONTINUE 

500 FORMAT C IX , F7 . 2 , F8 . 1 , F7 . 1 , 2F6 . 2 , 2F10 . 1 , F12 . 1 , F8 . 2 D 
DO 58 1-25,85,20 

XLOC-C I-5D/80. 0+1.0 
TMAR-TCR I T-TWAC I D 

WRITEC6, 500DXL0C, TWAC I D/TMPC, TWCCI D/TMPC, YYCTWAC ID, YAD/YLD1CID, 
& YYCTWCCO. YAD/YLD2CID, HACID* HCNV, HCCI V HCNV, QPPC ID. TMAR 

58 CONTINUE 

WRITEC6, 501DTHTHK.UMT, TB 

501 FORMAT C//' THTHK- ' , F6 . 4 , 4X , * UMT- ' , F4 . 1 , 4X , ' TB- ' .F5.1D 
WRITEC6.903D 

903 FORMAT C ' 1' D 

IFCTYPE . NE. 'P'D GO TO 1234 
22 CONTINUE 
21 CONTINUE 
20 CONTINUE 

1234 PRINT 4 ,' R-REPEAT W/ SAME MAT ' ' L N-NEW MAT ' ' L Q: ' 

READ' CAD ' , RES 
IFCRES . EQ. 'R'DTHEN 
GO TO 554 

ELSE IF CRES .EQ. 'N'DTHEN 
GO TO 555 
END IF 


1000 END 
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****** SUBPROGRAM 

SUBROUTINE DATAIN C EA , AA , P A , AK . YA . MET AL3 

REAL K. TMPC503 , EMPC5Q3 , EACO : 15 3 , AACO: 153, PACO: 153. AKCO: 153. 
& YACO: 153.AMPC503.PMPC503. AKPC503.YMPC503 
CHARACTER*50 METAL, NAME. JUNKC23, RES , FORMER . NAM 
FORMER- ’START' 

101 PRINT*.’ ENTER MATERIAL NAME: ' 

READ ' CA3 ’ .METAL 

1 CL0SEC93 
CLOSEC 103 
CLOSE 013 
CLOSE 023 
CLOSE 043 
CLOSEC153 

0PENC9, FILE-'MLIST’ , STATUS- ’ OLD' 3 
OPENC 10, FILE-' MUST’ , STATUS-' OLD’ 3 
OPENCU. FILE-' HOLD' . STATUS- ' OLD ' 3 
OPENC 12 , FILE-’ HOLD' . STATUS- ’ OLD ' 3 
0PENC14, FILE-' CHOICE’ 3 
OPENC 15, FILE-’ CHOICE' 3 
REWIND 9 

READC9, ' CA3’ 3CJUNKCI3, 1-1,13 

DO 2 1-1,100 

READC9, 700, END-33NAME 

WRITEC12, 701 3NAME 

IFCFORMER . EQ. ’ % ' 3THEN 

WRITEC15.7013NAME 

END I F 

FORMER-NAME 

700 FORMAT OX, A3 

701 FORMAT C2X.A3 

IF CNAME .EQ. METAL3 GO TO 5 

2 CONTINUE 

3 PRINT*,' LIST OF AVAILABLE MATERIALS’ 

REWIND 14 

READC14, ’ C A3 ’ 3CJUNKCJ3. 0-1, 13 

DO 30 1=1,100 

READ C 14 , 700 , END-31 3 NAM ' 

30 WRITEC* .7013NAM 

31 PRINT*,’ R-RETRY, N-ENTER NEW MATERIAL, Q-QUIT’ 

PRINT*,' FOLLOW YOUR SELECTION WITH A CARRAIGE RETURN ’ 
READ’ CA3 ’ , RES 

IF CRES .EQ. 'R' 3 THEN 
GO TO 101 

ELSE IF CRES .EQ. ’ Q' 3THEN 

STOP 

END IF 

WR I TEC 12 , 701 3METAL 

PRINT*,’ YOUNG' ’ S MODULUS SECTION’ 

6 PRINT*,' INPUT TEMP. C~1 TO EXIT3 ’ 

READ* ,TMPK 

IFCTMPK .LT. 03 GO TO 8 

PRINT*,’ INPUT YOUNG' ’S MODULUS ’ 

READ*, EL 

WRITE C 12 , 8003 TMPK.EL 
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800 FORMAT C2X , F6 . 1 . E2Q . 9) 

GO TO 6 

8 WRITE (12,801)' COEFFICI ENT OF THERMAL EXPANSION' 

801 FORMAT C2X,A) 

PRINT*,' COEFFICIENT OF THERMAL EXPANSION SECTION' 

9 PRINT*,' INPUT TEMP. C-l TO EX I TO ' 

READ* , TMPK 

IFCTMPK .LT. 0) GO TO 10 

PRINT*,' INPUT COEFFICIENT OF THERMAL EXPANSION * 
READ*, AL 

WRITE C 12 , 800) TMPK, AL 
GO TO 9 

10 WRITE (12,801)' POISSON RATIO' 

PRINT*,' POISSON RATIO SECTION' 

11 PRINT*.' INPUT TEMP. C-i TO EXIT) ' 

READ* .TMPK 

IFCTMPK .LT. 0) GO TO 12 
PRINT*,' INPUT POISSON RATIO ' 

READ*, PR 

WRITE (12,800) TMPK, PR 
GO TO 11 

12 WRITE (12, 801) 'THERMAL CONDUCTIVITY' 

PRINT*,' THERMAL CONDUCTIVITY SECTION' 

14 PRINT*,' INPUT TEMP. (-1 TO EXIT) ' 

READ*, TMPK 

IFCTMPK .LT. 0) GO TO 15 

PRINT*,' INPUT THERMAL CONDUCTIVITY K ' 

READ*, K 

WRITE (12,800) TMPK, K 
GO TO 14 

15 WRITE (12,801)' YIELD STRENGTH' 

PRINT*,’ YIELD STRENGTH SECTION' 

16 PRINT*,' INPUT TEMP. (-1 TO EXIT) ’ 

READ* , TMPK 

IFCTMPK .LT. 0) GO TO 17 
PRINT*,' INPUT YIELD STRENGTH ' 

READ*, YS 

WRITE (12,800) TMPK, YS 
GO TO 16 

17 WRITE ( 12 , 801 ) ' V 
REWIND 11 

READCli, ' (A)' XJUNKCI), 1-1,1) 

DO 7 1-1,100 
READCli, 700, END-1) NAME 
WR I TEC 10, 701 ) NAME 
7 CONTINUE 

PRINT*,' FILE LENGTH EXCEEDS DO LOOP LIMIT' 

GO TO 1 
5 CONTINUE 

DO 20 1-1,50 

RE ADC 9, 802 , ERR-21 )TMP( I ) , EMP( I ) 

20 CONTINUE 

802 FORMAT ( IX , F6 . 1 , E20 . 9) 

21 CONTINUE 

CALL TEMP (TMP , EMP . I -1 , EA) 



DO 22 1-1,50 

READC9 . 802 , ERR-23JTMPC I ) , AMPC I ) 

22 CONTINUE 

23 CALL TEMPCTMP, AMP. 1-1, AA) 

DO 24 1-1.50 

READC9, 802 . ERR-25DTMPC I ) , PMPC I ) 

24 CONTINUE 

25 CALL TEMPCTMP. PMP. 1-1, PAD 
DO 26 1-1,50 

READC9 , 802 , ERR-27DTMPC I 5 , AKPC I ) 

26 CONTINUE 

27 CALL TEMPCTMP, AKP, I-l.AK) 

DO 28 1-1,50 

READC9 , 802 , ERR-29JTMPC I ) , YMPC I ) 

28 CONTINUE 

29 CALL TEMPCTMP, YMP. I-l.YA) 

1000 END 


***• SUBPROGRAM *■*■*•• 

SUBROUTINE TEMPCT, P, M, CD 

....... F I NDS THE POLYNOMIAL COEFFICIENTS OF UP TO A THIRD 

....... DEGREE POLYNOMIAL FIT TO A MATERIAL PROPERTY VS, TEMP. 

....... USING METHOD OF LEAST SQUARES 

REAL TC50). PC 50), AC 10. 10), SUM1 CO : 15} , BCO: 15),CC0: 15) 

INTEGER CMAX 
COMMON/DF/DF 
DO 5 1-0,15 
SUM1CD-0.0 
BCD-0. 0 
5 CCD-0.0 

IFCM . EQ. DTHEN 
CC0)-PC1) 

RETURN 
END IF 

IFCM .LE. 3 DTHEN 

CMAX-M 

ELSE 

CMAX-3 

END IF 

DO 10 Kl-1 , CCMAX-1) *2 
DO 10 JJ-l.M 

SUM1 CK1)»SUM1CK1) + CTC JJD/DF) * * K1 
NOTE DIVIDING FACTOR * 

10 CONTINUE 

DO 20 K2-0 , CMAX-1 
DO 20 JJ-l.M 

BCK2)-BCK2)+CTCJJ)/DF)* ‘K2*PCJJ) 

DIVIDING FACTOR * 

20 CONTINUE 



DO 30 I-l.CMAX 
IM-I-1 

DO 30 J-l.CMAX 
ACI, JD-SUMiCIMD 
30 IM-IM+1 

ACI, 1D-FL0ATCMD 
NN-10 

CALL GAUSSC A , 3 , C , NN, CMAX D 

RETURN 

END 


***** SUBPROGRAM ******* 

SUBROUTINE GAUSS CA , B , X , ND , NO 

REAL ACND, NDD, BCND), XCNDD, AINVC5Q,50D 

DET-1.0 

DO 1 I - 1 , N 

DO 1 J-1,N 

IFCI .EQ. JD THEN 

AINVCI, 0*1.0 

ELSE 

AINVCI. JD-0.0 
ENDIF 

1 CONTINUE 

DO 9 I PASS-1, N 

I MX* I PASS 

DO 2 I ROW- I PASS, N 

IF CABSCACIROW, IPASSDD . GT . ABSCACIMX, IPASSDDD THEN 

I MX- I ROW 

ENDIF 

2 CONTINUE 

IFCIMX .NE. I PASS} THEN 
DO 3 I COL* 1 , N 
TEMP-AINVCI PASS , I COL} 

A I NV C I PASS , IC0L3=AINVCIMX, I COLD 
AINVCIMX, I COLD=TEMP 
IFCI COL .GE. I PASSDTHEN 
TEMP-AC IPASS, ICOLD 
ACIPASS, ICOL)-ACIMX, ICOL) 

AC I MX , ICOD-TEMP 
ENDIF 

3 CONTINUE 
ENDIF 

PIVOT-ACIPASS. I PASS) 

DET-DET * P I VOT 
IFCDET . EQ . O.ODTHEN 
WRITEC* , 10} 

STOP 

ENDIF 

DO 6 ICOL-l.N 

AINVC IPASS, ICOLD-AINVCI PASS , I COO/P I VOT 
IFCI COL .GE. I PASS} THEN 
ACIPASS, ICOLD - ACIPASS, ICOLD/PIVOT 
ENDIF 

6 CONTINUE 
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DO 8 I ROW-1, N 
IFCIROW . NE. IPASSD THEN 
FACTOR-ACIROW, IPASSD 
END IF 

DO 7 I COL-1 , N 

IF CIROW .NE. IPASSD THEN 

AINVC I ROW, ICOLD-AINVCI ROW , ICOLD-FACTOR'AINVC IPASS. I COLD 
ACIROW, ICOLD-ACIROW, ICOLD-FACTOR*ACIPASS. ICOLD 
END IF 

7 CONTINUE 

8 CONTINUE 

9 CONTINUE 
DO 4 I-l.N 
XCID-0.0 
DO 4 J-l , N 

4 XCID-XCID+AINVCI , JD’BCJD 
RETURN 

10 F0RMATC5X, ' ERROR IN GAUSS THE MATRIX IS SINGULAR',/, 

& 10X, ’ PROGRAM TERMINATED’ D 
END 


* * * * SUBPROGRAM * * * * * 

FUNCTION YYCT, CD 
REAL CCO: 15D 
COMMON/DF/DF 
YY-0.0 
DO 10 1-0,5 

10 YY-YY+CCID*CT/DFD**I 
END 


....... SUBPROGRAM ******* 

FUNCTION RADIUSCXD 

• FINDS RADIUS AS FUNCTION OF AXIAL LOCATION 

• DISTANCE MUST BE IN METERS 

XI-X/0.0254 

• CONVERSION TO INCHES 

IF CXI .GE. 1.7D THEN 
RADIUS-CO. 2253*XI+0. 08199D* .0254 
RETURN 

ELSE IFC XI .LE.- 0.60D THEN 

RAD I US-C-0 . 39825 * X I +0 . 07705D * 0 . 0254 

RETURN 

ELSE 

XI -X I +0 . 60 

RADIUS-0. 316- . 4266 *XI + . 73 *X! *XI-. 6132* XI * *3 + 

& .2828* XI * *4-.0273*XI * *5-. 0231 * XI * *6+ . 00587 * X I * *7 
RADIUS-RAD I US*. 0254 
RETURN 
END IF 
END 
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***** SUBPROGRAM ***•••• 

SUBROUTINE SOLVEC XA . XC , EPS , I MAX , ROOT , F) 

Xl-XA 
X3-XC 
I * * 0 

F1-FCX1) 

F3-FCX3) 

D-X3-X1 

1 X2-Xl-D*Fl/CF3-Fl) 

F2-FCX2) 

IFCD . LT. EPS) THEN 
R00T-X2 
I TER* I 
RETURN 

ELSE IF CABSCFCX2)) .LT. EPS)THEN 
R00T-X2 
I TER- I 
RETURN 

ELSE IF Cl .GT. IMAX) THEN 
WR I TEC * , 10) I , X2 , F2 
STOP 
END IF 

I FCF1 * F2 .LT. O.OTHEN 

X3-X2 

F3-F2 

FI- . 9 * FI 

ELSE IF CF2 * F3 . LT . O-OTHEN 

X1-X2 

F1-F2 

F3- . 9 * F3 

ELSE 

WR1 TEC * > * ) ' NO ROOT IN I NTERVAL ‘ , XI , X3 , ' IN STEP', I 
STOP 
END IF 
D-X3-X1 
1*1 + 1 
GOTO 1 

10 F0RMATC5X. 'NO ROOT HAS BEEN FOUND I N’ . 1 4 , ' ITERATIONS' , / , 
5X , ' THE LATEST VALUES ARE X- ',F10.5, 

' AND FCX)- ' , El 1 . 4) 

END 


• SUBPROGRAM ******* 

FUNCTION FMACMA) 

COMMON/AR/AR 
REAL MA 

FMA-AR-C1 .0+0.2 *MA*MA)* *3/Cl . 728 *MA) 
END 


* SUBPROGRAM ******* 

FUNCTION ROMBCF,A,B»KMX, EPS , IPRNT) 

REAL TCO : 15 , 0 : 15 ) , RCO s 15,0: 15) , DX, A , B, EPS , SUM 

INTEGER K,M,NPTS, KMX, IPRNT 

K-G 

DX-CB-A) 

TCQ,0)-0.5*DX*CFCA)-FCB)) 


1 


K-K+l 
NPTS-2* *K 
DX-DX/2 
SUM-0.0 

DO 2 I-l.NPTS-1.2 
SUM-SUM+FC A+ I * DX) 

2 CONTINUE 

TCK,0)-TCK-1.0)/2.0+DX*SUM 
DO 3 M-O.K-l 

TCK.M+l)-TCK.M)+CTCK.M)-TCK-i,M))/C4.0**CM+l)-l.G) 

IFCK .GE. 2 .AND. M .LI. K-2)THEN 

RCK.M)-CTCK-1.M)-TCK-2.M))/CTCK.M)-TCK-1,M))/4.0**CM+1) 

IF CABSCRCK,M)-1.0) .GT. 0.75 )THEN 
WRITEC* . 10)RCK.M).TCK-1,K-1) 

IFCIPRNT . EQ. 1)G0 TO 5 
ROMB-TCK-1 , K-l) 

RETURN 

ENDIF 

ENDIF 

3 CONTINUE 

IFCABSCTCK, K)-TCK-1 , K-l )) -LT. EPS)THEN 

romb-tck, K) 

IFCIPRNT .EQ. 1)G0 TO 5 
RETURN 

ELSE IFCK .LT. KMX)THEN 

GO TO 1 

ELSE 

WRITEC * . 1 1 )KMX , TCK, K) , TCK-1 , K-l ) 

ROMB-TCK, K) 

IFCIPRNT .EQ. 1)G0 TO 5 

RETURN 

ENDIF 

OUTPUT TABLE 

5 CONTINUE 

WRITEC*. 12X1. I-O.K) 

DO 6 I-O.K 

WRITEC*, 13)1, CTCI ■ M) , M-0, I ) 

6 CONTINUE 

WRITEC*. 14X1. I-O.K) 

DO 7 1-2, K 

WRITEC *.15)I,CRCI,M), M-0. 1-2) 

7 CONTINUE 
ROMB-TCK, K) 

RETURN 

FORMATS 

10 F0RMATCT20, ’ WARNING FROM ROMBERG INTEGRATION FUNCTION’,/, 

+ T20, ’ ROUND-OFF ERROR FLAG - *,F9.5./, 

+ T20 , ’ MOST RECENT DIAG. ROMBERG TERM- ’.E12.6,/, 

+ T20, ’CALCULATION HALTED AND THIS VALUE RETURNED’) 

11 F0RMATCT20, 'WARNING FROM ROMBERG INTEGRATION FUNCTION',/, 

+ T20, 'THROUGH ORDER ’,12,' THE REQUIRED ACCURACY’,/, 

+ T20, 'COULD NOT BE ACHIEVED. THE LATEST VALUE',/, 

+ T20, ' RETURNED TO THE CALLING PROGRAM WAS - '.E12.6,/, 
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+ T2G , ' AND THE PREVIOUS VALUE IN TABLE IS * ',E12*60 

12 FORMATC/// , TIO . 'THE ROMBERG TABLE’ ,//, T5 , ' LEVEL ->',/, 

+ IX, ' ORDER' ./, 1X.9C6X, 12, 6X)0 

13 FORMATC1X. 12, 1X.9CF13.9, 1X00 

14 FORMATC//, TIO, ' THE TABLE OF THE ROUND-OFF ERROR FLAGS’.//, 

+ TIO, 'LEVEL ->',/, IX, 'ORDER' ,/,4X, 15C3X, 12,3X00 

15 FORMATC1X, 12, 1X.15CF7.4, 1X00 
END 


* SUBPROGRAM ******* 

FUNCTION TOPCS0 

COMMON/FT/ST , TR , TW , TO , DRDZ , MA , DELMA . STEP , RAD , CF , DSDTH , T I NF , DELTW 
B-TO/TW 

C-CTO-TINFO/TW 

TOP-S* *7/Ci • O+B*S-C*S*S0 

END 


* SUBPROGRAM ******* 

FUNCTION BOTCS0 

COMMON/FT/ST, TR,TW, TO, DRDZ, MA, DELMA, STEP, RAD, CF, DSDTH, T INF, DELTW 
B-TO/TW 

C-CTO-TINFO/TW 

BOT-S* *7* Cl- O-S0/C1 • O+B* S-C*S*S0 
END 


SUBPROGRAM ******* 

FUNCTION TRAP CF,A,B,IDV0 

DX-CB-A0/IDV 

X-A 

TERM1-FCX0 
TRAP-O. O 
DO 10 I-O, IDV 
X-X+DX 
TERM2-FCX0 

TRAP-TRAP+G. 5*DX* CTERM1+TERM20 
TERM1-TERM2 
10 CONTINUE 
END 


** SUBPROGRAM ******* 

FUNCTION FTHTACTHTA0 

COMMON/FT/ST, TR,TW, TO, DRDZ.MA, DELMA, STEP, RAD, CF, DSDTH, T INF, DELTW 
REAL MA 

FTHTA-CF/2 . 0* Cl . 0+ DRDZ * DRDZ 0 * *0 . 5-THTA * 

& CC2. O-MA* MA+DSDTHO/ CMA+O . 2 *MA* * 30 * DELMA/ STEP+ DRDZ /RAD 0 
END 


** SUBPROGRAM •***••* 

FUNCTION F1PHICPHI0 

COMMON/FT/ST, TR, TW, TO, DRDZ, MA, DELMA , STEP , RAD , CF , DSDTH , T I NF , DELTW 
REAL MA 

F1PHI— PHI “Cl • O-MA * MAO/ CMA + O. 2 *MA* * 3 0 * DELMA/STEP + DRDZ/RAD 
& -DELTW/STEP/CTO-TWO 
END 



** SUBPROGRAM * 

FUNCTION F2PHICPHO 

COMMON/ FT/ST , TR , TW , TO . DRDZ , MA , DELMA . STEP , RAD , CF , DSDTH, T I NF . DELTW 
REAL MA 

F2PHI »ST * C CTR-TWD/CTO-TVO ) * C 1 • + DRDZ * DRDZ ) * * . 5 
& -PHI*CCl-0-MA*MA5/CMA+0.2*MA* * 3} 'DELMA/STEP+DRDZ/RAD 
& -DELTW/STEP/CTO-TVOD 
END 
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Temperature distribution in the nozzle wall 


We can write the general expression for heat transfer 
through the nozzle liner wall: 


This can be 




rewritten 


r 

I k (T) dT 
T 

w,a 

ln < r 0 / r i> 
2 ft 4kL 

as 


(El) 


q k In (r/r x ) = -k 21T*L (T - T^ a ) 


Now solve for r. 


let 


then 


T - T = — — 

w , a 2 Th k A.L In 

~3k 

2 fhk In r. 

T = K In (t) + T w _ a 


(E2) 


Thus, the temperature through the nozzle liner wall has 
been shown to be a logarithmic function of the radial dis- 
tance . 
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APPENDIX F 


Comparison of Thermal and Pressure Stresses 

The highest stresses, and therefore the lowest factors 
of safety, occur at the throat. It is in this region that 
we should compare the pressure stresses with the thermal 
stresses. For the inside surface of the nozzle liner we 
have thermal stress defined by Equation 10 in the main text. 
The pressure stress is defined by Equation 14. 

At throat conditions for the material TZM with a liner 
thickness of 1.5 mm we have the following values: 
oc = 0.554 X 10“ 5 m/m K 
E = 0.23 X 10 12 Pa 
T = 400 Kelvin 

r^ = 0.00569 m 

r = 0.00719 m 
o 

= 0.3 

P = 7.29 X 10 6 Pa 

P = 13.8 X 10 6 Pa 
c 

Using Equations 10 and 14: 


( . ) . . .= 3.93 X 10 

v 0 ,i' thermal 


8 


t 


and 




6,i'pressure 


= 2.10 X 10 
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The ratio of thermal stress to pressure stress is 18.7. 
We conclude that the nozzle design will be based largely on 
thermal stress considerations. 


APPENDIX G 


Pre-blowdown Nozzle Failure 

Before the blowdown process is begun the coolant side 
of the nozzle liner will be subjected to the full 2000 psi 
coolant pressure, while the air side will experience a hard 
vacuum. Two modes of failure must be considered in this 
case: compressive failure and buckling failure. 

Compressive failure. For this analysis we use the same 
triaxial state of stress treatment used for the operational 
condition except that we remove thermal stresses and 
stresses due to air pressure. In this way we can calculate 
a Factor of Safety. The results of this analysis are shown 
in Figure FI. It can be seen that the minimum Factor of 
Safety is 4.18 at the location furthest upstream of the 
throat. Compressive failure, then, in the pre-blowdown 
condition is not a problem. 

Buckling failure. For this analysis we employ buckling 
failure theory for a right cylindrical shell. A more accu- 
rate model would be to use analysis for a conical shell. 


However, a cylindrica 

1 shell 

is 

less resistant 

to failure 

than 

a conical shell, 

and so 

we 

will obtain 

a 

very conser 

tive 

result. If this 

result 

is 

acceptable 

for 

a cylinder 

we can be quite sure 

it will 

be 

acceptable 

for 

a cone. 
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The equation used for this analysis comes from Faupel 

[ 26 ] : 


> . = E 

Crlt 4 R 3 (1 - > 2 ) 


(FI) 


where R and t are the local nozzle radius and thickness, 
respectively, and X is Poisson's ratio. The results of 
this analysis are shown in Figure F2 . The minimum critical 
pressure for buckling is seen to be 6300 psi. Since 2000 
psi is the maximum pressure the liner wall will be subjected 
to, buckling failure is not a problem. 
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Figure FI. Compressive stress for pre-blowdown condition 
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